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Abstract. The calculated effects of interstitial hydrogen on the elastic properties of 
a-iron from our earlier work are used to describe the H interactions with homogeneous 
strain fields using ab initio methods. In particular we calculate the H solublility in Fe 
subject to hydrostatic, uniaxial, and shear strain. For comparison, these interactions 
are parametrised successfully using a simple model with parameters entirely derived 
from ab initio methods. The results are used to predict the solubility of H in spatially- 
varying elastic strain fields, representative of realistic dislocations outside their core. 
We find a strong directional dependence of the H-dislocation interaction, leading to 
strong attraction of H by the axial strain components of edge dislocations and by 
screw dislocations oriented along the critical (111) slip direction. We further find a 
H concentration enhancement around dislocation cores, consistent with experimental 
observations. 



PACS numbers: 62.20.-x, 61.72.Bb, 64.30.Ef, 81.40.-z 
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1. Introduction 

The influx of liydrogen into many alloys and steels, either during manufacture or service, 
is unavoidable and leads to undesirable consequences. One of these is the lowering 
of the failure stress, leading to embrittlement and fracture at unpredictable loading 
conditions [H [2]. Some proposed mechanisms of H-embrittlement of iron include the 
hydrogen-enhanced decohesion (HEDE) mechanism [Sj H] in which H is postulated to 
cause a loss of cohesion between the metal atoms leading to failure at interfaces, H- 
vacancy effects O [6], |7] in which vacancies containing H may order themselves along 
critical slip directions, and hydrogen-enhanced localised plasticity (HELP) [HI El 
where the onset of plasticity with loading occurs at a lowered stress as a result of the 
H shielding of repulsive interactions between dislocations [9]. Due to the increased H 
concentration at dislocations [H], this shielding effectively reduces the spacing between 
dislocations [lOj. Support for the HELP hypothesis in explaining the phenomenon 
of dislocation coalescence, and ultimately, crack advancement at reduced loads, is 
provided by the experimentally-observed increase of dislocation mobility [10]. Various 
models of H-embrittlement are discussed in the context of recent experimental results 
by Kirchheim [12]. H is very mobile in Fe, with a diffusivity of 10~^ — 10~''m^/s as 
determined from a wide range of experiments ^3], and by calculations [14j. Therefore, 
it can thus be assumed to be present in all microstructural regions. In real materials 
with defects, H is most prevalent at dislocations in the case of single crystal Fe, and at 
grain boundaries in polycrystalline Fe [15] . 

The above theories assert that H degrades the material, but the details of the 
mechanisms are still under debate [161 |T7]. Nevertheless, H-dislocation interactions have 
long been implicated in the modification (whether degradation or enhancement) of the 
strength properties of iron [161 113 and recent experiments [18] confirm the importance 
of increased H-concentration near dislocations and other low-energy trap sites. 

Compounding the confusion caused by the existence of multiple theories of H- 
embrittlement is that these various mechanisms have mainly been proposed and 
discussed separately and rarely together [191 H] as would be desirable for an unbiased 
assessment of their relative importance under various conditions. The advantage of 
theroretical approaches is that the effects can be investigated independently. However, 
a major challenge for achieving a theoretical understanding of these mechanisms lies 
in the vastly different length scales involved: the large number of atoms required as a 
result of the slow spatial decay of dislocation strain fields and small H concentrations 
versus the highly-localised modification of the electronic structure around H. As a result, 
simulated systems need to be extremely large in order to account for the long-ranged 
elastic distortions induced by dislocations. Taketomi et al [20] conducted a molecular- 
statics study with empirical potentials of the H distribution around an edge dislocation in 
Fe. Ref. pi] describes a kinetic Monte Carlo study using empirical potentials to examine 
diffusion and trapping of H in a bulk Fe crystal containing screw dipoles. Clouet et al [22] 
have used empirical potentials combined with elasticity theory to model the interaction 
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of C atoms with edge and screw dislocations in Fe. However, empirical potentials have 
limited and largely untested transferability to these systems [23], and more accurate 
studies using first-principles techniques are required. The current drawback with first- 
principles methods is that only small system sizes can be treated, leading to idealised 
systems or conditions. The static calculation of the energy barrier for H-diffusion in 
bulk Fe [m 123], the study of the stress distribution around H in Fe pi], or a first- 
principles molecular-dynamics study of H-diffusion in a very small (16 atom) Fe unit 
cell [25], are typical of the types of problems that can currently be treated using first- 
principles methods. Therefore, first-principles density functional theory is frequently 
combined with other techniques in order to examine effects requiring both atomistic 
accuracy and a large number of atoms. Examples of recent studies include a study of the 
microstructural evolution of tungsten irradiated with He [26], parametrising elasticity 
theory for studying solute-dislocation interaction [27], and making extrapolations to 
continuum-scale properties such as changes in the shear modulus, due to the effect of H 
on the electronic density of states of Fe [28] . 

The problem of interstitials interacting with dislocations has long been recognized 
to be important and has been studied analytically [H] |29l |30l [31l [321 [33l |33] within 
continuum theory. A review of early work on the subject can be found in Ref. [35] . 
These studies recognized that the orientation of the dislocation with respect to the unit- 
cell axes can have dramatic effects on the solubility and henceforth on the interstitial 
solute concentration around dislocation cores. However, while they are able to capture 
the effects of edge dislocations, whose strain fields are dominated by axial stresses, by 
not going beyond hnear elasticity, the models with hydrostatic expansion due to the 
solute fail to yield any interaction between solutes and screw dislocations, whose purely 
shear-stress fields do not interact with the hydrostatic component of the strain field 
of the interstitial solute and studies thus ascribe the interaction of interstitials with 
screw dislocations purely to the tetragonal distortion which they introduce. An earlier 
study [M] did obtain interactions within linear elasticity theory with screw dislocations 
but only because their model assumes a tetragonal distortion field generated by the 
interstitial solute. Very recently, the interaction of C with dislocations was investigated 
by ab initio methods [36] in a model within linear elasticity theory which accounted for 
tetragonal distortions. This limitations of the models in properly describing trapping 
of interstitials by screw dislocations has long been noted and ascribed to a non-linear 
elastic, or "modulus effect" in the hterature [371 133 EH] but studies of solute-dislocation 
interaction have mostly tended to ignore this effect. A recent exception is the work in 
Ref. [39] which studies the solute-dislocation interaction at the continuum scale. 

The goal of this study is to use ab initio results from our previous study [IQ] in order 
to make connections with the proposed theories of H-embrittlement for a-Fe. Our results 
of the influence of interstitial H on the elastic properties ^40J are used to parametrise 
an elasticity-theory model for the behaviour of H within strain fields representative 
of dislocations and to predict concentration enhancements around typical dislocation 
strain fields. 
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In Sec. [2] we describe the details of our DFT calculations used to calculate the 
solubility of H under various strain fields. The elasticity-theory model of H solubility 
under spatially-homogeneous strain is introduced in Sec. [3] and the predictions are 
compared with the explicitly-calculated results. The model is further apphed in Sec. H] 
to study the solubility in strain fields oriented in different directions with respect 
to the unit-cell axes in order to establish a correspondence with locally-modified H 
concentrations in the vicinity of edge and screw dislocations. We conclude with a 
discussion in Sec. [5l 

2. First-Principles Calculations 

The details of the calculations are the same as in our previous study [IQ] and are 
therefore only summarised here. Our spin-polarised first-principles density functional 
theory calculations were performed with the VASP [HI US], HS] code. We used the 
projector augmented-wave method [HJ |15] and considered the 3p electrons of Fe 
as valence electrons. The generalized-gradient approximation (GGA) in the PW91 
parametrisation [46j was used for the exchange-correlation functional. A plane-wave 
basis with a cutoff of 500 eV and supercell geometry were used. A F-centred k-point 
grid equivalent to 18x18x18 for the two-atom basis bcc unit cell was used for the 
Brillouin-zone sampling except for the 128-atom Fe cell where it was the equivalent of 
20x20x20. 

We focus on tetrahedral interstitial H because (i) at zero stress we find it to 
be 0.13 eV more stable than the octahedral position, in agreement with previous 
DFT studies [TU [2l], and (ii) we expect it to dominate the mechanical properties at 
ambient temperatures as there are only half as many octahedral interstitial sites. In 
all calculations, the cells with H were taken as cubic in the unstrained state, which is 
a good approximation, for low concentrations of H, for the determination of the elastic 
parameters PU]. The tetragonal distortion, which is usually accounted for explicitly is 
studies of octahedral interstitials such as C in a-Fe [SU [36], is smaller in the case of 
interstitials occupying tetrahedral sites. In addition, the cubic approximation makes 
sense as the occupancy of interstitial sites by H is expected to be random. 

For the results presented here, we used different H concentrations which we arrived 
at using different supercells (Table II in Ref. |1Q]). For some of the cases studied 
(Sec. I3.3.2P it was necessary to calculate the Fe and Fe-H cells at the same stress. Equal 
stresses in the two systems were maintained by enforcing lattice parameters such that 
the stresses in the two systems are equal according to Eq. [TH which is valid for small 
strains. The first-principles determination of stress, from the DFT program, served as 
a further check that the systems had the same stress. 
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3. Strain dependence of H solubility 

3.1. Solution energy 

In this section, we describe the effect of different types of strain on H solubility. Our 
goal is to examine the effect of H on typical dislocation strain fields and conversely, 
the latter's impact on H solubility. Instead of simulating dislocations explicitly in 
computationally very demanding large supercells, we have examined the solubility of 
H in different types of homogeneous strain fields which we later use to estimate the H 
solubility and henceforth the concentration profile in realistic dislocation strain fields. 

Specifically, we have calculated the solution energy for tetrahedral interstitial H 
in supercells subject to hydrostatic, or uniaxial or shear strain. These applied strain 
tensors mimic the strain components of different dislocations as will be discussed in 
SecH 

To arrive at the solution energy W ofnB. atoms, at each value of applied strain, we 
subtracted from the total energy of the cell with N Fe atoms and n H atoms, Epe^Hn, 
the total energy Ep^j^ of a pure iron cell with Fe atoms, and half the energy of a free 
H2 molecule multiplied by n: 

W = Epe^H^ - Epe^ - n^EH2- (1) 

The solution energy can be interpreted in different ways depending on the strain 
or stress state of the pure Fe reference system with respect to the system with H. We 
have defined three regimes on which to concentrate as the strain of the Fe-H system 
is varied: (i) the pure Fe volume matching that of the Fe-H system at each strain, (ii) 
the two systems having the same stress tensor at each strain, or (iii) having the same 
strain tensors. Cases (i) and (iii) would be described by the change in Helmholtz free 
energy with constant volume and temperature, and constant strain and temperature 
respectively, with H, while case (ii) would be described by the change in the Gibbs free 
energy where temperature and pressure are maintained steady. Case (i) would reflect 
the energy to be gained or lost by inserting |if2 into a strained Fe lattice without further 
changing the average volume per Fe atom. It would be the expected boundary condition 
for H in an perfect infinite Fe lattice. It would also correspond to regions where H is 
present in already-expanded regions of the metal, e.g. near a grain boundary or a 
dislocation. Cases (ii) and (iii) would involve other types of confinement, for instance, 
in the case of H locally present at crack tips, the applied boundary conditions would 
determine whether any eventual failure is stress- or strain-controlled. The boundary 
conditions chosen are important for a variety of systems, affecting, for example, the 
microstructure of alloys |17] or the activation energy of dislocation nucleation 

3.2. Parametrisation of the H solution energy 

Using our DFT-obtained material parameters (Table [1]), we have parametrised the 
solution energies (Eq. [T]) as functions of strain and H-concentration. We label all 
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quantities pertaining to the Fe-H or the pure Fe system with an x or superscript 
respectively: for example, V"^ for the equilibrium volume where Y = x for the Fe-H 
system and Y = for pure Fe. For a system Y the total energy at an arbitrary reference 
volume V^'^-^ corresponding to a stress tensor cr^'^^ and infinitesimal strain tensor e with 
respect to V^'^-^ is (c./. Eq. 6 in Ref. [ID]) comprised of three parts: a reference term, a 
stress term, and a strain term, 

(y-^ e) = 0) + v^^f Yl + ^ E crMy''^)^^^^M. (2) 

ij ijkl 

The first term is the total energy at the reference volume V^'^-^ , the second term is the 
stress part which vanishes for V^'^^ = ^ and the last term is the elastic deformation 
energy about V^^^ . The expression ([2]) equivalently arises from an expansion of the 
elastic energy from up to third order in strain defined with respect to . Explicitly, 
the reference term in Eq. [2] is equal to 



2 ^ ' dV 



2 



= E^iV\0) + ^{V-f-V^)\ (3) 

where -B^ is the bulk modulus of system Y at its equlibrium volume V"^ . The second 
term in Eq. O assuming a hydrostatic stress corresponding to a pressure P^'^-^ required 
to reach the reference volume, is 



ij ij 

= - P^'^^V'^^Y^n 
i 

= £J^(V"/-V-)5^e„, (4) 

i 

where P^^-^ = —B^ (y^^-^ — V^^ jV^ was used. This relation for P"^"^^ makes use of 
linear elasticity, i.e. that doesn't vary, in Eq. HI which is a reasonable assumption 
for low concentrations. In the third term of Eq. |2l the elastic parameters at the volume 
V^^^ are equal to 



V=VY 

(5) 



'ijkiy.' )^ yY 

' V=VY 

where the volume derivative has been rewritten in terms of the hydrostatic strain 
r]{V) = {V-V^) IV" . 

Despite their generality, Eqs. [3}|5] are difficult to work with for the Fe-H system, 
due to the complex strain-dependency of the Cj^^i- We prefer therefore to express all 
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quantities in terms of small variations from x = and V^°, and we take the reference 
volume as = and expand about as = + xdV/dx (for small x). 
With these expansions, we arrive at a general form of Eq. |2]for the Fe-H system that we 
will use throughout the applications presented in the following sections. The first term 
of Eq. [2] (z.e. Eq. [3]) is thus equal to 



1 



dE 



dx 



d^E 



d^E 



1 o d^E 



dE 



VO,x=0 

dV d'^E 



+ 



dx"^ 



dxdV 



VO,x=0 

dV 



VO,x=0 



X' 



dx 



(6) 

, d^E 
dx"^ ' 



dx ' dx dxdV ' 2V^ 
where here and in what follows, the limits of evaluation of the derivatives are omitted 
for clarity. The first derivative with respect to volume vanishes as the evaluation limits 
correspond to pure Fe. The second term of Eq. [2] {i.e. Eq. H]) also vanishes because 
is the system's equilibrium volume. For the third term, the elastic parameters (Eq. [5]) 
can be written as 



where the derivative represents the variation with respect to concentration about 



X 



(7) 



and with strain included implicitly. These derivatives were determined in Ref. 
are given in Table [H 



0] and 



Table 1. Parameters used in the solubility expressions [TOllTS] and in subsequent 
sections. All derivatives are evaluated at V'^ and a; = 0. 



Quantity 


value 


A E 


0.22 eV 




194.2 GPa 


A V 


4.5 A3 


V^N 


11.352 A3 


dB°/dr] 


-1075 GPa 


dC^Jdr] 


-1492 GPa 


dC^^/dr] 


-669 GPa 


dClJdri 


-761 GPa 


dB*°ydx 


-153 GPa 


dClf/dx 


-316 GPa 


dC{f/dx 


-72 GPa 


dClf/dx 


-141 GPa 



After defining the relative volumes, stresses, and/or strains, of the Fe-H and Fe 
systems corresponding to each of the cases outlined above, we use Eqs. [T]and[2]to derive 
strain-dependent solution energies for different relative strain and/or volume conditions 
for the Fe and Fe-H systems. 
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3.3. Applications of solution energy model to limiting cases 

3.3.1. Same volumes The reference volume of the Fe-H system, , is taken as the 
reference volume of the pure Fe system, resulting in a contribution from the stress term 
(Eq. H]). The pure Fe system, at a reference volume and strain e^, applying the 
equations of the previous section (Eq. [2}|lj), is thus described by a total energy 



RO / rl\/\ ^ T/^ rl\/ . - T/^ - - 



2V« \ dx I V" dx 

i ijkl 



kl 

(8) 



where the elastic parameters are (Eq. E}, 



V° dr] 



The solution energy (Eq. [T]) for a cubic system is then given by the contributions to 
Eq. [2] where we neglect all terms of order except for those in Eq. [6l which are required 
for explaining the discrepancies at zero strain between the various concentrations: 



+ 



1 
2 



2N (9x2 Q^gy 
dCjf _ dC^ 
N dx df] 



Av] fe" 2 + 2 + 

I \ XX ' yy ' zz ] 

dx drj J yy xx^zz ^ ^yy^zz) 

( V^dC^ _ ^AV^^ fe" 2 ^ ' + e" 

^ c/x rfr^ V ^ ^ ^ ^ 



(10) 



In the above, we have used the definition of concentration as the ratio of the number 
of H atoms, ra, to the number of Fe atoms N, x = n/N, to rewrite the concentration 
derivative of V in terms of the volume expansion AV per H and A^: 4^ = 4^4^ = Al^A^ 

^ ^ ax an ax 

and similarly for the concentration derivative of E, ^ = AEN where AE = j^ — \Eh2- 
The ratio /N is the (constant) volume per Fe atom at equilibrium. In Eq. [TOl the 
x-dependent terms - 0{x'^) in - are given only to highlight the higher-order effects 
at zero strain; otherwise the expression constitutes a universal formula valid for all low 
concentrations. We obtained a value of A£'=0.22 eV by calculating the \H2 solution 
energy in the largest unit cell (128 atoms) studied. The expression for the solution energy 
can be simplified for particular cases of applied strain. For the case of hydrostatic strain 
(- = exx = eyy = e^z, Eq. HO] becomes 

W'^y'^yn = AE + ^^ + xAV^ - SB'AVe 
2N ox^ oxoV 



+ 1 _rr AV t\ (11) 

2\N dx dr] ' ^ ^ 
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Here, we have not associated the quantity | ^ — — with because the slopes of the 
energy-strain coefficients as a function of pressure are not equal to those of the stress- 
strain coefficients [m |50] used to calculate the bulk modulus. However, the numerical 
values of the respective slopes are not very different [101 • The simplified, linear in 
strain, form of Eq. [TTl W/n = AE — SB^AVe, is a well-known result of the literature 
of continuum mechanics [TTl \51\ . 

The corresponding expressions for uniaxial strain e = exx and shear strain e = exy = 
dyx are 

' 2N dx^ dxdV \N dx dr] J ^ ^ 

A combination of hydrostatic/uniaxial/shear strain corresponds to combining 
expressions [TTIfl3l as can be seen from Eq. [TOl 

By evaluating the coefficients of Eq. [TOl apart from those of the x-dependent terms 
which are discussed below, using the values given in Tab. HI we obtain a parametrisation 
of the H solution energy for small arbitrary strain, pressure, and H-concentration. To 
that end, we performed DFT simulations with the Fe and Fe-H systems at the same 
volume and varied the amount and type of strain. We subtracted the two sets of energies 
and the energy of (Eq. [T]). The parametrised forms (Eqs. [TTIfT^ . using the DFT- 
calculated parameters listed in Table [H accurately describe the DFT-obtained solution 
energies, as can be seen from Figure [H The parametrised curves in all of the figures 
were generated by assuming the dilute limit [lima;_j.o (^)] of Eq- HOl such that the 
contribution from the x-dependent parts of Eqs. [TOlfTSl is zero. 

The variations of the calculated solution energies at zero strain among the various 
concentrations can be accounted for by evaluating the terms a^ff + ^ g^gy + 
{c.f. Eq. [6] and the normalised form in Eq. [TOl) explicitly at each concentration as 
summarised in Tab. HI The cross-term was calculated for x=1.8% and this result 
(-0.557 eV/A) was used throughout. Using a graph similar to Fig. [1^, but with the 
horizontal axis expressed in terms of the strain e° of pure Fe, the sum of the ffist two 
terms of Eq. [TO]was predicted using the zero-strain data (equal to a volume of from 
there. This can be done since W{x)\^o=q = AE + 2^^i-- These sums are given as the 
second column of Table [2j The total sums including the cross-term, stated in the fourth 
column, agree very well with the DFT data given in the last column. The cross-term was 
also calculated using data from the 5.9% concentration and it was a bit smaller (-0.426 
eV/A) than the one used for Table |2l leading to a contribution xAV g^^y =-0.11 A eV 
and thus a total E^{V^, 0)=0.141 eV, for the entry 'sum' at x = 5.9%, thus accounting 
for the minor deviations at the higher concentrations between the two last columns of 
the table. The deviation of the intercept from the zero-concentration limit of 0.22 eV 
for all the graphs shown in Fig. [T] can thus be explained. 



Ab initio parametrised model of strain- dependent solubility of H in a-iron 



10 




Figure 1. [color online] DFT-calculated solution energies along with parametrised 
form for (a) hydrostatic strain, (b) uniaxial strain, (c) shear strain, for different H- 
concentrations. All strains are component (length) strains. 
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Table 2. Energy terms from Eq. [3] contributing to the intercept at zero strain for the 
case of equal volumes (Fig.[T]). The first two energy columns are from expression [B] and 
were calculated using DFT as described in the text. Their sum is under the heading 
sum. The last column is the DFT data. 



x{%) 




(eV) 




sum 


DFT W\^^^^ieY) 


0.8 


0.246 




-0.023 


0.223 


0.229 


1.8 


0.238 




-0.048 


0.190 


0.192 


2.7 


0.249 




-0.066 


0.183 


0.190 


3.6 


0.251 




-0.089 


0.162 


0.175 


4.0 


0.250 




-0.093 


0.157 


0.173 


5.3 


0.278 




-0.123 


0.155 


0.177 


5.9 


0.255 




-0.149 


0.106 


0.136 



The solution energy as a function of strain and concentration can also be derived 
by a systematic expansion of the total energy in terms of strain and concentration |27j . 
In that work, the authors presented a systematic expansion of the total energy in 
terms of strain and concentration. The resulting derivative of pressure with respect to 
concentration yields a contribution to the solution energy which is linear in hydrostatic 
strain, in agreement with our findings, and with the earlier work in Refs. pT| |5T] which 
assumed a continuum model throughout. Using the notation of Ref. [27] to express 
our concentration and strain-dependent elastic parameters requires considering up to 
the third derivatives of the total energy with respect to concentration and strain but 
otherwise yields equivalent results. 



3.3.2. Same stresses or same strains As discussed at the end of Sec. 13. H another 
possibility of evaluating the solution energy (Eq. [1]) is to compare the Fe-H and pure Fe 
reference systems for the same stress or strain tensors. Therefore, we expand the total 
energy of each system about its equilibrium volume. 

Using Eq. [21 the same stress an, for the two systems corresponds to 



(Tr, 



nk^k 



(14) 



in the Voigt notation for m,n,k from 1 to 6. In the case of hydrostatic strain e^j = dije, 
this leads to e° = e^. An expansion of the stress (energy) to second (third)-order in 
strain leads to contributions in the solution energy of powers of which we ignore. If 
instead of equal stresses, the systems have equal strain tensors, the relation is simply 



'ij- 



We obtain, function of Fe-H hydrostatic strain e 



6ije and shear strain 



-xy 



-yx 



6, the following expression for the solution energy per H: 

h xAV 1 

2N dx"^ dxdV 2 



W/n = AE 



9 



B^AV± 



dxdV 
iV dx 



xAV^ 



S\ (15) 




Figure 2. [color online] DFT-calculated solution energies per H atom for various 
concentrations along with parametrised forms (Eq. llSp as a function of shear strain 
where the reference pure Fe system has (a) the same strains, and (b) the same stresses. 
The strain on the horizontal axis is for the Fe-H system. 



The upper (lower) signs are for equal strains (stresses). Eq. [T5]is sufficiently general 
such that expressions can be found for mixed cases, e.g. equal hydrostatic stresses 
but equal shear strains. The expression [15] is compared with directly explicit DFT 
calculations in Figure [2] for e = and varying 5 for the cases of (a): equal strains, and 
(b): equal stresses. The higher curvature for x=2.7% in (b) is due to the C44 data (Fig. 
4 in Ref. |10]) lying below the line of best fit, thereby yielding a larger contribution from 
the total derivative dClf /dx. 

As with the previous case of equal volumes, the intercept at zero strain in Fig. [2] 
can be accounted for by explicitly calculating the zero-strain terms appearing in Eq. [151 
along with an extra term (third column of Table [3]). The comparison with the DFT 
data, in Table |3l is very good for small concentrations and deviates for the larger one, 
by a similar amount as the case of constant volumes. As in that case, discussed at the 
end of Sec. 13.3. H the deviation can be attributed to variation of the cross-term with 
concentration. 
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Table 3. Energy terms contributing to the intercept at zero strain (Fig. [2]) for the 
case of equal strains or stresses. The first energy column is the same as that in Table [2] 
while there is now an additional term (second energy column in this table) appearing. 
The sum of the first two columns is under the third column, sum, which should be 
compared with the DFT data in the last column. 



x{%) 


Table [2] sum 




sum 


DFT VF(x)|,.^o(eV) 


2.7 


0.183 


0.030 


0.213 


0.220 


4.0 


0.157 


0.045 


0.202 


0.211 


5.9 


0.106 


0.068 


0.174 


0.202 



4. H solubilities in realistic dislocation strain fields 
4.I. Anisotropy of H solubility 

The parametrisation of the H solution energy in terms of strain presented in the previous 
section enables us to determine the H solubility in typical dislocation strain fields. We 
will focus on the particular case of a screw dislocation lying along (111), which is the 
predominant orientation of a screw dislocation in a bcc metal [52]. Within isotropic 
dislocation theory [3H], the strain field of a screw dislocation is described in cylindrical 
(r, 6, z) coordinates by its only nonzero component as 

^ez = ^ze = -y— (16) 
47rr 

where 6 = 2.45 A is a Burgers vector equal to the length of |(111) for a dislocation 
lying along (111). 

The solution energy for pure shear strain at the same volume (Eq. [T3|) and at the 
same stress/strain (Eq. ITSj) is independent of the choice of shear direction as long as it 
is oriented along the cube axes i.e. if e%y = e^^ = 6, = e^^ = 6, or e^^ = e^y = 6. A 
screw dislocation line oriented along (001) would result in shear components oriented 
along the cube axes, specifically with e^^ = ef^ 7^ 0, e^^ = ^ty 0; e^y = e^^. = 0. 
Referring to Fig. IH we see that H is attracted towards the case of the screw dislocation 
in the case where the pure Fe reference system has the same shear strain as the Fe-H 
system (Fig. [2^), but repelled in the case of same stresses (Fig. |2]d) and same volume 
(Fig. [ID . 

However, such an attractive behaviour can vary significantly with the direction 
of applied shear strain [53]. Therefore, the solubility equations [TTIIT3] and [15] were 
transformed to rotated coordinates corresponding to orientations other than the (100)- 
type considered in previous sections. By rotating the entire coordinate system, 
{x,y,z) {x',y',z'), other directions z' for the dislocation line are achieved. The 
transformation of strains from the reference (unprimed) to rotated (primed) coordinate 
system can be performed using Euler angles (see e.g. Ref. [38]): e' = TeT""" where T 
is the transformation matrix for transforming the reference coordinates x into the new 
frame x', such that x' = Tx. We expressed the strains e^'z' = ^z'x' in the primed system 
in terms of the unprimed system and inserted these values into the general expression 
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Figure 3. [color online] Variation of solution energy for 2% shear strain ey'z' (see 
inset) in the case of equal volumes for a hypothetical screw dislocation line/Burgers 
vector oriented along different z' directions (labels and {hkl) in inset). The horizontal 
axis is the variation of the solution energy as a function of the angle (j) which represents 
the choice of x' and y' axes within the plane normal to the Burgers vector, (p = has 
been arbitrarily chosen to correspond to the x' = {khO) direction (see inset). The 
dotted line is the unstrained value of the H solution energy. 



for the solution energy Eq. [10] for the case of equal volumes and its equivalent general 
form for the cases of equal strains/stresses. 

For the case of purely hydrostatic strain, the solution energy is independent of 
the coordinate system chosen because the trace of the strain tensor is invariant under 
any unitary transformation. Figure [3] shows the variation of the solution energy in the 
case of equal volumes, for ey'z'=2% where Eq. [10] has been rotated onto other axes. 
These axes are shown in the inset where the dislocation is taken to lie along z'. The 
z' direction may be associated with the Miller indices hkl of a plane. 0, the angle that 
x' makes with the (khO) direction, also defines the direction of y'. In Fig. [3] only the 
^y'z' — ^t'y' components are applied. For the full description of the screw dislocation 
strain field eg/^', the e^,^, = e^,^, components are also required. In the case of the (111) 
screw dislocation line, the solution energy as a function of the full strain field e^/^' 
remains fiat, independent of shear direction 0. In contrast with Ref. [53], we do not find 
asymmetry of our solution energy W about zero strain as a function of applied shear 
strain, but we cannot provide any insight into this discrepancy because the details of 
the model used in Ref. [53] are lacking. 

The curves in Fig. [3] while corresponding to the same magnitude of shear, do 
not necessarily correspond to the same distance r away from the dislocation core, as 
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Figure 4. [color online] Variation of H solubility energy as a function of distance 
from the (111) screw dislocation core for the cases of equal stresses, equal volume, or 
equal strains for the Fe-H and pure Fe systems. The dotted line denotes the zero-strain 
solubility. The strain field of the screw dislocation was represented by Eq. [1^1 



the Burgers vector [c.f. Eq. [TB]) is different for dislocation lines lying along different 
directions. The radial dependence of the H solubility was calculated for the case of a 
(111) dislocation for the three cases described in SecO The solubility corresponding to 
the complete strain tensor Eq. [16] was calculated from Eq. [10] in the case of equal volumes 
or analogous equations for the cases of equal stresses and equal strain. Regarding 
the case of equal stresses, we have not modelled the solubility of H in the stress 
field of a dislocation per se but its solubility at individual stresses at each point r, 
meaning that each infinitesimal annular element is not in mechanical equilibrium with 
its surroundings. However, this is a first approximation in a self-consistent treatment of 
H solubility - see e.g. Ref. ^1]. The results, shown in Fig. [H show that H is attracted 
to the core for equal shear strains or equal volumes and repelled for equal shear stresses. 

4.2. Concentration of H around dislocation strain fields 

Motivated by the strong strain-solubility effects, we calculated the local modification 
of the H concentration by strain fields for the case of equal volumes. We consider the 
change in free energy F as a result of adding H to Fe under strain. The free energy is 

F = Ef,.h-TS, (17) 

where T is the temperature, and S is the total entropy of the Fe-H system. In order 
to isolate the effect of adding H to the already-strained Fe lattice, we assume that 
the configurational entropy of the H is independent of that of the Fe lattice: i.e. 
S = Sh + Spe- Epe-H is the total energy of the (strained) Fe-H system. It can be 
rewritten in terms of the solution energy W {c.f. Eq. [1]) as Epe-H = Epe + W — fin, 
where we have set the energy per H as — /x + ^Eh^, in contrast to Eq. [T] where it was 
set to \Eh^. 
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The free energy is minimized with respect to n: 

SF _^_dW r^9S _^ 
5n dn dn 

We first focus on the case of a hydrostatic strain field. For simphcity we truncate 
Equation [TT] from the previous section to hnear order in strain: 

W = W{e) = n{AE - 3B°AVe). (19) 

The expression for the entropy can take on various forms, in the presence of various 
types of H-H interactions, all equivalent at low concentrations. In what follows, n 
is the number of H, N is the number of matrix (Fe) atoms, M = 6N is the total 
number of tetrahedral interstitial lattice sites, and c is the occupation concentration 
defined as c = n/M. The simplest form for the entropy, 5*0, arises from an unrestricted 
concentration of H, corresponding not to a lattice but to an ideal gas. In this case, the 
entropy derivative is ^ = — fcln-p. The concentration yielded by Eq. [TS]for the case of 
hydrostatic strain up to linear order in the strain takes the Maxwell-Boltzmann form of 

c = coexp [3B°AVe/kT] (20) 

where Cq = exp [{—AE + fi)/kT]. This reproduces the result from the Cottrell picture of 
enhanced concentration around a dislocation [TT] . A more realistic form of the entropy 
takes into account the presence of a lattice and enforces a maximum occupancy of one H 
per tetrahedral site. In this case the entropy is 5*1 = —k [nlnfj + (M — n)ln (l — -p)] 
and the entropy derivative is ^ = —kin The concentration using this form follows 
Fermi-Dirac statistics: 

c = {exp [(A^ - - W^AVe) /kT] + l}"^ (21) 

However, the maximum occupancy of one H per tetrahedral site is not realistic as the 
H-H interaction in Fe is repulsive: We have calculated the H-H interaction as a function 
of distance and find strong H-H repulsion up to third-nearest neighbours where it was 
0.13 eV. For a realistic representation of the statistics, up to and including third nearest- 
neighbour H-H spacings should be strongly unfavourable in our model. In the extreme 
case, these sites are left unoccupied, corresponding to infinite repulsion. The expression 
for the entropy derivative for tetrahedral interstitials in a bcc lattice with total exclusion 
of first nearest-neighbours, both first and second-nearest-neighbours, and up to third- 
nearest neighbours is given in Ref. jES] as 

dSinn _ 6n _ z-^'^~ 3M)(5n2 - 36nM + 72M^) 

dn ~ ~ ^zin + MzY ^~ (n - QMfM ' 

dS2nn , , n{M - n) 
— km 



dn (M - 4n)2 

dS2,nn 



dn 



khi 



n{M - ?,n){M - 2nf^ 

_ {M-An){M-hnY \ ^ ^ 
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respectively. In the limit of low concentration, all of the expressions for the entropy (5*0, 
Si, Sinn, S2nn, S^nn) are the Same, the differences being of O(c^). Similarly, at small 
strains, all of the solutions for c are identical. Therefore, the H-H interaction can be 
ignored for predicting H concentrations around dislocations. 

A numerical example of the extent of the modification of the H-concentration in 
bcc Fe by applied hydrostatic strain is shown in Figure |5^. In this low-concentration 
regime, all forms of the entropy are identical. In this example, the background (zero- 
strain) atomic H concentration (n/N) was taken to be equal to 3x10^^, in line with 
typical densities of H in iron [2J. At each temperature studied, /i was adjusted in 
order to achieve this background concentration. From the figure, it is clear that the H- 
concentration can be dramatically increased by orders of magnitude, even at very small 
strains. The general trend of increased H concentration under tensile strain supports 
the interpretation of H solubility as a volume effect. The increase in concentration is 
more pronounced at lower temperatures, in agreement with the low thermal desoprtion 
temperatures, ca 400 K, indicated for H trapped at dislocations |18]. We take our results 
of Fig. |5^ to be representative of edge dislocations which have a significant hydrostatic 
strain component [38]. With shear strain applied in the 6'z' direction where z' is the 
(111) direction, and regardless of the direction of 6', we obtain the curves in Fig. |5]d. 
The effect is less pronounced than for hydrostatic strain, because the solubility energy 
is quadratic in the shear strain (c./. Eq. [131 whose coefficients change upon rotation to 
this coordinate system, but whose functional form remains the same). 



H cone (H/Fe) H cone (H/Fe) 




hydrostatic strain shear strain 



Figure 5. [color online] Concentration of H as a function of (a) hydrostatic strain, 
and (b) shear strain egi^' where the dislocation line z' is taken as the (111) direction, 
for a chemical potential giving a background concentration of 3 x 10~^ H per Fe atom 
at each temperature. The case of constant volumes for Fe and Fe-H was assumed. 

By combining component strains in our solubility expressions, we are able to predict 
H concentration profiles in the full strain fields of edge, as well as screw dislocations. 
H concentration profiles are shown for two common slip systems in a-Fe in Fig. [61 It 
is expected that symmetry reduction will result from the consideration of tetragonal 
distortions, as for the case of C in Fe in Ref. [3^, and by the inclusion of dislocation 
anisotropy [Ml |56] . 
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Figure 6. [color online] Concentration contours of H around an (a) edge dislocation, 
with the tensile region at y < and (b) screw dislocation, at 300 K, for a background 
concentration of 3 x 10~^ H per Fe atom. The case of constant volumes for Fe and Fe-H 
was assumed. The shaded circle denotes a region of radius three times the Burgers 
vector at the core. 

5. Conclusions 

We have derived a model for the H point-defect - dislocation interaction (solubility) in 
a-Fe and parametrised it with ab initio- derived quantities from a previous study |I0] . 
Although, for simplicity, our model assumes a hydrostatic lattice expansion due to H, 
we find not only an effect on the solubility from axial strain, but from shear strain 
as well. In order to get an interaction with screw dislocations, it was necessary to 
go beyond the H- induced size effect and consider the modulus effect [371 [38l [39] . 
From our models we find, in agreement with existing experimental studies [151 IlH], 
and a recent atomistic study [36], that typical strain fields found near dislocations can 
trap substantial amounts of H. We further show that the local H-concentration may be 
enhanced by several orders of magnitude in hydrostatic strain fields, observed in edge 
dislocations, and this enhancement identifies a realistic regime for the consideration 
of H concentrations in the atomic % range as in our previous study jlO]. In the case 
of shear strain, which constitutes the strain field of screw dislocations, despite our 
simplifications of hydrostatic expansion of the lattice by H, we also find attraction of H 
if the screw dislocation is oriented along the (111) direction. However, the concentration 
enhancement for shear strain is not as great as in the case of hydrostatic strain. We 
find, in general, that the areas of enriched local concentration tend to occur mainly at 
low temperatures, in the hundreds of K. Similar conclusions regarding the enrichment 
of C in a-iron as a function of strain and temperature were reached in a recent atomistic 
study [36]- Screw dislocations along (111) are the predominant types of dislocations in 
bcc metals driving plasticity J52j and the concentration enhancement which we find near 
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screw dislocations is consistent with various theories and observations of H trapping at 
dislocations. Our findings support the strong role played by dislocations [151 IIH] in 
describing the effect of H on the mechanical properties bcc Fe. 
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